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We present an accurate equation of state for water based on a simple microscopic Hamiltonian, with only 
four parameters that are well-constrained by bulk experimental data. With one additional parameter for the 
range of interaction, this model yields a computationally efficient free-energy functional for inhomogeneous 
water which captures short-ranged correlations, cavitation energies and, with suitable long-range corrections, 
the non-linear dielectric response of water, making it an excellent candidate for studies of mesoscale water 
and for use in ab initio solvation methods. 



I. INTRODUCTION 

The emergence of several macrosco pic phases of wa- 
ter with distinct microscopic structure*^ from relatively 
simple molecular interactions places this liquid at the 
forefront of interesting unsolved problems in the study 
of condensed matter. The structure of water around 
microscopic objects differs significantly from the bulk, 3 
and such effects play a critical role in the structure of 
proteins^ and in chemical reactions at catalyst surfaces.^ 

Current computational approaches to systems sensitive 
to liquid structure most often employ molecular dynam- 
ics simulations. Ab initio molecular dynamics,^ which 
treats all the valence electrons in the system quantum 
mechanically, is relatively accurate but prohibitively ex- 
pensive for all but the smallest of systems. Hybrid ap- 
proaches that combine electronic structure methods for 
part of the system with classical molecular dynamics sim- 
ulations for the fluid can handle larger systems, but re- 
quire empirical models for both the electron-fluid cou- 
pling and the classical pair potentials for the fluid. These 
molecular dynamics methods are inherently expensive 
due to the sampling required for thermodynamic aver- 
ages which, when coupled to an electronic system, ne- 
cessitates repeated electronic structure calculations. In 
addition, the thermodynamic integration required to cal- 
culate free energies, which are necessary for analyzing 
chemical reaction pathways, significantly exacerbates the 
cost of such methods. 

Efficient theories for the equilibrium properties of liq- 
uids, on the other hand, deal directly with average den- 
sities instead of individual molecular configurations. In- 
tegral equation theories give accurate structures of in- 
homogeneous fluids,^ but still prove relatively expensive, 
particularly for estimating free energies. The most direct 
approach to free energies is classical density-functional 
theory, a method based on approximations to the ex- 
act free-energy functional of the liquid^ which has the 
added advantage of being readily coupled to electronic 
density-functional theory within the framework of joint 
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density-functional theoryP The most accurate, currently 
available functionals for polar molecular fluids such as 
water P^HSl however, rely on direct correlations (from 
neutron scattering or computer simulation) at each tem- 
perature and pressure of interest, restricting their effi- 
ciency and applicability. 

This work addresses the need for a computationally 
efficacious microscopic theory of water that is capable 
of providing accurate free energies under inhomogeneous 
conditions without the dependence on fluid structure 
data. The strategy is to identify a simple effective mi- 
croscopic Hamiltonian which (a) reproduces the equation 
of state for homogeneous water and (b) is readily repre- 
sented by a free-energy functional even in the inhomoge- 
neous case. 

Statistical Associating Fluid Theory,^ based on 
Wertheim's thermodynamic pertubation theory,^ is one 
such approach which has been successfully applied to the 
study of vapor-liquid interfaces,^ with model parame- 
ters for water determined from the equation of stated 
However, the predictions of this theory for the inhomo- 
geneous fluid have not yet included quantities of interest 
in solvation methods such as pair correlations, cavitation 
energies and dielectric response, partly due to the relative 
complexity of the model Hamiltonian. Below, we develop 
an alternate simpler Hamiltonian based upon microscopic 
intuition about hydrogen bonding, and we demonstrate 
that the resulting functional (also based on Wertheim 
theory) leads to a relatively accurate free-energy descrip- 
tion of inhomogeneous water, especially given the sim- 
plicity of the underlying model. 



II. MODEL MOLECULAR HAMILTONIAN AND THE 
EQUATION OF STATE FOR WATER 

Within the constraints of condition (b) above, a nat- 
ural starting point would be the standard approach of 
perturbation about the hard- sphere fluid, for which Fun- 
damental Measure Theory-SuEIl provides a highly accu- 
rate functional. The hard-sphere diameter required to 
reproduce bulk properties can be inferred from the ex- 
cluded volume in the equation of state, and fits^ to ex- 
perimental data suggest a value that strongly decreases 
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FIG. 1. (a) Tangentially bonded hard-sphere model for liquid 
water: O (oxygen) sphere with two V spheres in contact dia- 
metrically opposite to the H (hydrogen) sites (b) Residual for 
the equation of state fit to experimental datap^^ com- 
pared to the semi-empirical Jeffery- Austin equation of stateP^ 



with temperature and is ~ 3.3A at 298 K. This is clearly 
incompatible with the almost temperature-independent 
~ 2.8A location of the first peak in the experimentally 
observed oxygen-oxygen radial distribution! 2 ^ 

This incompatibility stems from the discrepancy be- 
tween the close-packed coordination of the hard-sphere 
fluid and the tetrahedral coordination favored by wa- 
ter. Water prefers the formation of open tetrahedral net- 
works at lower temperatures, which leads to empty space, 
"voids", within cages of water molecules, as manifested 
by the temperature-dependent excess excluded volume in 
the equation of state. Consequently, we propose a refer- 
ence fluid consisting of a compound object (FIG. [lja)): 
a hard sphere of radius Ro at the O (oxygen) site with 
smaller spheres of radii Ry at two void sites V placed in 
contact (at a distance aov — Ro + Rv) along two of the 
conjugate tetrahedral directions to the hydrogen bond di- 
rections. For our model, we take the O-H distance to be 
1 A with a tetrahedral H-O-H angle, as in the frequently 
employed SPC/E interatomic potential modelP^The ge- 
ometry of this compound object is chosen to encourage 
closest approach along the hydrogen bond directions. 

Our ansatz for the intermolecular Hamiltonian is the 
repulsive pair potential corresponding to the tangentially 
bonded hard-sphere trimer of FIG.Qa), perturbed by an 
isotropic attractive pair potential U a {r) between the O 
sites. The equation of state of this fluid is well approxi- 
mated by 
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where the first three terms correspond to the bonded 
hard-sphere equation of stated and the final term is the 
mean-field contribution from the as yet undetermined at- 
tractive perturbation U a (r), with k = — J drATTr 2 U a (r). 

The bonded hard sphere equation of state is based on 
Wertheim perturbation theory^ about the hard sphere 
mixture, consisting of density n of 0-spheres and 2n 
of T^-spheres. The pressure of this reference system is 



p HS = 3nT + pgg, where we separate and collect the 
0(n) ideal gas parts to elucidate the connection with the 
density functional (Ml) . For p^ s , we employ the accurate 
generalization^ of the Carnahan-Starling excess pressure 
to hard sphere mixtures 
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2Ry)n and n 3 — -g-v-^o 
fundamental measures. 



3n, 7ii = (-Rq + 2R v )n, n 2 = ^(R 2 



2Ry)n are the uniform fluid 



First order Wertheim perturbation for the bonding 
constraints accounts for the fixed O-V separation and 
not the V-O-V angle; nonetheless it has been shown to 
well approximate the equation of state of objects with 
this geometry.^ We accumulate its contribution at 0(n) 
into the first term of (Jl]): this exactly corrects the ideal 
gas mixture value of 3nT to the rigid-molecule ideal gas 
value of pid = nT. (We use this fact later to restore the 
bond angle constraints in the intramolecular geometry of 
the inhomogeneous fluid.) The remaining contribution 
of this perturbation, the third term of (Jl]), corrects the 
excluded volume effects of the hard sphere mixture to 
account for the O-V distance constraints. There, 
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is the contact value of the O-V radial distribution in the 
hard sphere mixture with i?hm = RoRv l&ov- 

As motivated earlier, the temperature dependence of 
the exclusion volume is a critical feature of the equation 
of state for water .^Because the location of the first peak 
in the O-O radial distribution does not change apprecia- 
bly with temperature, we attribute this dependence to 
changes in the radii of the V spheres, modeled as a de- 
creasing function Rv{T) — Ry\~ T ^ Tv to qualitatively 
capture the effect of the empty spaces in the open tetra- 
hedral network. This leads to a model equation of state 
|l| with only four adjustable parameters (Ry\ TV, k, 
and Ro), which we fit to experimental data for the bulk 
liquid and vapoi^ including data for the supercooled 
liquidP 

The root mean-square error in the ratio of the pressure 
to the ideal-gas pressure, p/nT, is 4.8 x 10~ 2 for the cur- 
rent 4-parameter fit, which compares very favorably with 
2.9 x 10 -2 for the standard semi-empirical Jeffery- Austin 
equation of state^ (comparison in FIG.JlJb)), especially 
considering the fact that the latter fit employs more than 
twice as many (~ 9) adjustable parameters. Beyond pro- 
viding a reasonable fit to the equation of state, the key 
advantage of the present work is that these results stem 
directly from a model microscopic Hamiltonian, which 
we exploit below to construct a theory for the inhomoge- 
neous fluid. 
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FIG. 2. Metropolis sampling of a canonical ensemble of pa- 
rameters, shown in all six projections of the four-dimensional 
parameter space. One hundred random samples (+) are 
drawn from the full set (•) for error estimation of all sub- 
sequent results; x's mark the optimum parameter set. 



To ensure that our model parameters are indeed inde- 
pendent and physically meaningful, we employ Bayesian 
error estimation following Ref. [301 EU Specifically, we 
generate a canonical ensemble of parameter sets (FIG. [2]) 
with a Metropolis walk in parameter space, where the 
residual is the 'energy' and the 'temperature' is 2Cq/N p 
where Co is the minimum residual and N p = 4 is the 
number of fit parameters. The optimum parameters with 
standard deviations thus estimated are 
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(1.290 ± 0.049)A 
T v = (258.7 ± 12.3)^ 

k = (1.805 ± 0.074) x 10 a KA 3 
Ro = (1-419 ±0.010)A. 



(4) 



The modest eccentricities of the ensemble slices in 
FIG. [2] indicate that the covariances of these parameters 
are nominal, suggesting that the parameters have physi- 
cal meaning rather than merely controlling a flexible fit 
function for the equation of state. 

An advantage of this ensemble-of-models approach, 
which we exploit below, is that one can estimate how 
well the fit to bulk data constrains all subsequent pre- 
dictions, including those for the inhomogeneous fluid, by 
evaluating those predictions for a sampling of the ensem- 
ble of parameters (indicated in FIG.[2]as +'s), instead of 
just the one optimum parameter set. 



III. MODEL FOR INHOMOGENEOUS LIQUID 

Capturing the behavior of the inhomogeneous fluid re- 
quires information beyond merely the integrated strength 
k of the pair-potential interaction U a (r). This work 
demonstrates that the simplest next step, including in- 
formation about the range of the interaction, suffices to 
capture surprisingly well the main features of the short- 
range correlations in the liquid. To this end, we employ 
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FIG. 3. Partial radial distributions from the density func- 
tional |6]| (bundle of thin red lines for the ensemble, with 
results for optimum parameters highlighted by dashed line) 
and Monte Carlo simulations for the model Hamiltonian (blue 
circles), compared to experimental pair correlations of water 
from Soper et alP^ (green dot-dashed line) 



the attractive-part of the Lennard- Jones potential 



Ua(r) 



-9k 



Mr) 1 



r > 2 1 / 6 a u 
r < 2 1 / e cr u 



1/4, 

(5) 

which has the correct long range r 6 tail for the 
orientation-averaged interaction of a dipolar fluid. We 
fit the range o\j to reproduce the bulk surface tension at 
298 K (based on calculations with the free-energy func- 
tional below), finding ajj = 2.62 A. We thereby introduce 
only one additional fit parameter in going to the inhomo- 
geneous fluid. 

To evaluate the viability of this simple model Hamil- 
tonian for describing the inhomogeneous fluid, we com- 
pute its pair correlation functions (for each of the state 
points for which experimental correlations were measured 
in Soper et alPS) directly with canonical-ensemble Monte 
Carlo simulations of 2048 molecules. The comparison 
(FIG. [3]) between the behavior of this model microscopic 
Hamiltonian (circles) and the actual experimental cor- 
relations in physical water (green dash-dot line) is re- 
markable given the highly simplified form for the model. 
Although the secondary peaks in the O-O correlations of 
this model Hamiltonian do appear more at the character- 
istic distances for a hard-sphere rather than at those for a 
tctrahcdrally bonded fluid, the temperature and pressure 
dependence of the locations and heights of the first O-O 
peak compare reasonably to water. Similarly, although 
the first two peaks of the O-H and H-H correlations are 
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fused into a single broader peak, the general location and 
particle content of these peaks are in reasonable agree- 
ment for such a simple model. These details could be 
corrected in future work by fitting perturbation pair po- 
tentials of zero integral AU a p(r) with a, /3 £ {0,H} to 
the experimental correlation data, but the focus of the 
present work is the quality of predictions which can be 
made from a simple microscopic model with very few 
adjustable parameters (five) constrained purely by the 
macroscopic data. 

Having established a short-ranged microscopic model 
Hamiltonian which reproduces relatively well the exper- 
imental correlations in water, we turn next to develop- 
ment of a corresponding free-energy functional. The form 
of this functional, 



$M = $id + $hs + ®b + - I n (U a * n ) , (6) 



mirrors the equation of state (JlJ, and is composed of 
the ideal gas free energy, hard sphere excess functional, 
bonding correction and mean field perturbation. 

We start with the exact grand free energy functional 
for the ideal gas of rigid molecules and thereby restore 
exact treatment of the intramolecular bond-angle con- 
straints; this approach is consistent with Wertheim the- 
ory since the latter yields the exact rigid-molecular ideal 
gas pressure pid = nT at 0(n) in the uniform limit. The 
free energy of the inhomogeneous ideal gas with chemical 
potential u in external site potentials V a (r) is written as 

$ id [V,n[V]] / dfn a (r)(V a (r) - i; a (r)) 

-(H + T) ( dfno(?), (7) 



employing ideal gas effective potential d 32 * 33 * ip a (r) for a £ 
{0,H,V} as the sole independent variables.^ Here, the 
site densities are dependent variables computed using 
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where ui £ 50(3)/Z2, where uio denotes the correspond- 
ing rotation for a vector, and where R a i are the site co- 
ordinates for a molecule in the reference orientation cen- 
tered at the origin with i = 1 for a = O and i £ {1,2} for 
a £ {H, V}. Note that we have simplified the above ex- 
pression using the Z 2 rotation symmetry of the molecule 
about its dipole axis. 

To treat the hard sphere mixture excess free energy 
^gg, we use the 'White-Bear mark II' version of funda- 
mental measure theory^ (incorporating Tarazona's ten- 
sor modifications^) 
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in terms of the scalar, vector and tensor weighted densi- 
ties ni = wf * no + wf * ny for i £ {0, 1, 2, 3, vl, v2, m2}, 
where 



f fn \ 1 1 »3(2-w 3 )+2(l-rt 3 ) log(l-rt 3 ) 

2n 3 -3n 3 +2n 3 '+2(l-n 3 ) 2 log(l-ra 3 ) 



and f 3 (n 3 ) = 1 - 



(See the comprehensive review by R. RotfPfor details.) 
Note that this functional corresponds exactly to the hard 
sphere excess pressure ^ in the uniform fluid limit. 

Next, $f, accounts for the tangential bonding con- 
straints on the hard-sphere exclusion effects; note that 
the contribution from Wertheim pertubation to the ideal 
gas part has been absorbed into the exact rigid-molecule 
ideal gas free energy <3>id. The Helmholtz-energy density 
for this term in the uniform fluid limit is determined from 
the third term of u\ to be — 2nTlog gov(°~ov), which we 
generalize to the inhomogeneous version 
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with the vector correction factor £ = 1 — \n V 2\ 2 / n 2 . We 
include this factor here following the spirit of Yu et al.J^ 
where £ was introduced in analogy with the occurence 
of the vector weighted densities in the hard sphere mix- 
ture functional in order to improve agreement with Monte 
Carlo calculations. Finally, the last term in (|6| describes 
the attractive perturbation potential within a mean-field 
picture. 

The partial radial distribution functions implied by the 
free energy functional as evaluated from its analytic 
second variational derivatives using the Ornstein-Zernike 
relation, are in excellent agreement with the Monte Carlo 
simulations (circles and corresponding curve in FIG.[3j). 
(The minor artifacts in the interior of the hard cores are 
caused by the bonding correction, whose inhomogeneous 
generalization is not perfect.) The small spread in these 
results with variation of parameters in the ensemble ex- 
emplifies how tightly the bulk data indeed constrain these 
predictions within the assumed model. 



IV. PREDICTIONS FOR THE INHOMOGENEOUS 
LIQUID 

To evaluate the predictions of the above density func- 
tional for the inhomogeneous fluid, we perform direct 
minimization of ^ using the nonlinear conjugate gradi- 
ents algorithm^ with the values of ipaif) on a discretized 
grid as the independent variables. The orientation inte- 
grals involved in evaluating the site densities from the 
site potentials ^ are discretized using quadratures on 
50(3)/Z 2 . The calculations presented below are per- 
formed on radial or planar d = 1 dimensional grids^ 
where the azimuthal symmetry simplifies the orientation 
quadrature from 50(3) = § 2 x §1 to § 2 , which we tesse- 
late using a recursively subdivided icosahedron! 3 ^ 
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FIG. 4. Energy of the vapor-liquid interface as a function of 
temperature, compared to the experimental values for surface 
tension.^ 
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FIG. 5. Comparison of density functional predictions with 
SPC/E molecular dynamics results^ for (a) Radial distribu- 
tion around spheres of radii 2, 4, 6, 8, and 10A that exclude 
the oxygen sites and (b) the variation of solvation energy of 
such spheres per surface area with radius. 



We find remarkable agreement with available data for 
the behavior and free energies of inliomogeneous aqueous 
systems, especially given that only bulk data, including 
surface tension at a single temperature, were employed 
in determining the limited number of parameters in the 
functional. For example, FIG. [4] compares our prediction 
of the temperature dependence of the interfacial energies 
with experimental data. Over the entire range of accessi- 
ble temperatures at ambient pressure, we find the experi- 
ment to lie within the relatively narrow variations within 
our ensemble of models. Moving beyond planar inter- 
faces, FIG.[5]explores the radial distribution around hard 



FIG. 6. Nonlinear dielectric response: variation of relative di- 
electric constant e with externally applied field Eo at ambient 
conditions in comparison to SPC/Epl and BJtP 3 molecular 
dynamics simulations. By construction, the linear dielectric 
constant matches experiment due to the long-range correc- 
tions. 



spheres and the variation of free energy of hard-sphere 
insertion with radius, and demonstrates that the predic- 
tions of our model are in qualitative agreement with the 
SPC/E molecular dynamics results.^ The contact densi- 
ties and the free energies from our model are somewhat 
higher than those from the SPC/E model results, a situ- 
ation which could potentially be improved in future work 
by including additional perturbation pair-potentials. 

In addition to the bulk and short-ranged correlations 
described above, a successful theory of solvation requires 
accurate dielectric response. Following Lischner et al.p^ 
we add a scaled mean-field long range electrostatic cor- 
rection 



IK 



A e (T) 



E 



Z a Zf3 / n a K * rip, 



(11) 



where the site charges Z a are taken to be the SPC/E 
valued and K = G 2( 1+ (q/ g y-j with G c = 0.33 is the 
Coulomb kernel with a high frequency (short range in 
space) cutofPl The prefactor A e (T) = 1 - T/(7.35 x 
10 3 K) serves to correct for dipole correlations beyond 
mean field, and is fit to reproduce the bulk dielectric 
constant at small field. Fig. [6] shows that the nonlinear 
response at high fields (which is not fit) is well captured 
by the interplay between <E> e and <£>id- 

Conclusion — We have constructed a computationally 
tractable free-energy functional for studies of inhomoge- 
neous water based upon a microscopic Hamiltonian con- 
strained by experimental data for the bulk equation of 
state. Following this approach gives a remarkably high- 
quality fit to the equation of state with only four tightly 
constrained parameters. With one additional parameter, 
the range of the model interaction, the resulting func- 
tional captures the free energies associated with inhomo- 
geneous systems such as the liquid-vapor interface and 
the embedding free energy of microscopic objects, as well 
as essential features of the partial radial distributions and 
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density profiles around microscopic objects. With long- 
range corrections, the model gives an accurate descrip- 
tion of the non-linear dielectric response. The model thus 
shows good promise for capturing the key quanties which 
require description in solvation studies. In future work, 
further details may be captured with suitable perturba- 
tion of the pair potentials constituting the underlying 
microscopic Hamiltonian. 
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